\section{Code for problem 3}
\label{code3}
\small
\begin{verbatim}
# The \mu function with the parameters given in the assignment
mu <- function(t)
{
  return(10^(-3)+10^(-4)*1.1^(30+t))
}

# The number of repetitions
N <- 1000
# A vector for the simulated V(0)'s
Vs <- c()
for(j in 1:N)
{
  # The number of steps we take to find a single repitition
  NrSteps <- 1000
  # The length of the steps
  h <- 10/NrSteps

  # A vector for the simulated interest rate
  r <- rep(0,1001)
  # The initial interest rate
  r[1] <- 0.044
  for(i in 1:NrSteps)
  {
    # We simulate the interest rate in the next step by using Euler's method
    Z <- rnorm(1)
    r[i+1] <- r[i]+0.2*(0.05-r[i])*h+0.1*sqrt(h)*Z
  }

  # A vector for the simulated reserve
  V <- rep(0,1001)
  # The reserve at time T is 0
  V[NrSteps+1] <- 0
  for(i in NrSteps:1)
  {
    # We find the reserve in the previous step by using Euler's method
    V[i] <- V[i+1]-((r[i+1]+mu(i*h))*V[i+1]-5+mu(i*h)*1000)*h
  }

  # We collect the reserves at time 0
  Vs <- c(Vs,V[1])
}

# we find the estimate for the reserve by taking mean of our simulated reserves at time 0
sum(Vs)/N

# We plot a histogram of the reserves at time 0
hist(Vs, main="Histogram of simulated V(0)'s", xlab="")
\end{verbatim}

\subsection{Alternative solution}
\label{code3b}
\small
\begin{verbatim}
# Integrated mortality
Lambda <- function(t)
{
  return(10^(-3)*t+10^(-4)*((1.1^t-1)/log(1.1))*1.1^30)
}

# Function to solve for death time
f <- function(t, u)
{
  return(-log(u)-Lambda(t))
}

# Calculates the discount factor for a given time.
calcDiscounting <- function(time, interest, stepSize)
{
  d <- 1
  for(i in 1:length(interest))
  {
    d <- d + d * interest[i] * stepSize

    if(i*stepSize > time)
      return(d)
  }
}

# Calculates the premium contribution of a death time
calcPremiumContrib <- function(deathTime, interest, stepSize)
{
  n <- round(deathTime/stepSize)
  s <- 0
  for(i in n:0)
  {
    s <- s - (s * r[i+1] - 5) * h
  }
  
  return(s)
}

# The number of steps to use
NrSteps <- 1000
h <- 10/NrSteps

# A vector for the calculated reserves
reserves <- c()

# The number of interest rates to simulate
N <- 100

# The number of death times used to estimate the reserve for each simulated interest rate
n <- 10000

for(i in 1:N)
{
  # Simulate the interest rate
  r <- rep(0,NrSteps+1)
  r[1] <- 0.044
  for(i in 1:NrSteps)
  {
    Z <- rnorm(1)
    r[i+1] <- r[i]+0.2*(0.05-r[i])*h+0.1*sqrt(h)*Z
  }

  # Holds the total reserve
  totalReserve <- 0
  
  # Calculates the value of the reserve if the insured survives
  valueOfSurvivel <- calcPremiumContrib(10, r, h)

  for(i in 1:n)
  {
    # Simulates the time of death by inverse transform method
    timeOfDeath <- uniroot(f, c(0, 1000), u=runif(1))$root
  
    # If the insured survived to time 10 the 'valueOfSurvivel' is added to the 'totalReserve'
    if(timeOfDeath >= 10)
    {
      totalReserve <- totalReserve + valueOfSurvivel
    }
    # Else the value of the death time is calculated
    else
    {
      totalReserve <- totalReserve + 
                      calcPremiumContrib(timeOfDeath, r, h) - 
                      1000/calcDiscounting(timeOfDeath, r, h)
    }
  }

  reserves <- c(reserves, totalReserve/n)
}

# The estimate and histogram are calculated
mean(reserves)
hist(reserves)
\end{verbatim}

\newpage 
